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Abstract. The object of this paper is a one-dimensional generalized porous media equation 
(PDE) with possibly discontinuous coefficient /?, which is well-posed as an evolution prob- 
lem in L (K). In some recent papers of Blanchard et alia and Barbu et alia, the solution 
was represented by the solution of a non-linear stochastic differential equation in law if the 
initial condition is a bounded integrable function. We first extend this result, at least when 
P is continuous and the initial condition is only integrable with some supplementary techni- 
cal assumption. The main purpose of the article consists in introducing and implementing a 
stochastic particle algorithm to approach the solution to (PDE) which also fits in the case when 
ft is possibly irregular, to predict some long-time behavior of the solution and in comparing 
with some recent numerical deterministic techniques. 
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1. Introduction 

The main aim of this work is to construct and implement a probabilistic algorithm 
which will allow us to approximate solutions of a porous media type equation with 
monotone irregular coefficient. Indeed, we are interested in the parabolic problem 
below: 

f d t u(t,x) = \d 2 x j{u{t,x)), te [o,+oo[, 

u(0, x) = uo(dx), i£l, 

in the sense of distributions, where uq is an initial probability measure. If uq has a 
density, we will still denote it by the same letter. We look for a solution of (11.11 ) with 
time evolution in X'(R). We formulate the following assumption: 
Assumption(A) 

such that /3|k + is monotone. 

(ii) /3(0) = and ft continuous at zero. 

(hi) We assume the existence of A > such that (/3 + Xid)(R + ) = (R+), id(x) = x. 
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A monotone function (3q : R — >■ R can be completed into a graph by setting /3n(x) = 
\fio(x-), (3q(x+)]. An odd function (3q : R — >■ R such that /3|r + = /3o|r+ produces in 
this way a maximal monotone graph. 

In this introduction, however (3 and (3q will be considered single-valued for the sake 
of simplicity. We leave more precise formulations (as in Proposition 12. 1 1 and Theorem 
12.81 1 for the body of the article. 

We remark that if (3 fulfills Assumption(A), then the odd symmetrized /3o fulfills the 
more natural 

Assumption(A') 

(i) /?o : R — >■ R is monotone. 

(ii) /?o(0) = and (3q continuous at zero. 

(iii) We assume the existence of A > such that (/3n + AicQ(R) = (R), id(x) = x. 
We define $ : R ->• R+, setting 




if u ^ 0, 



S(u) = < (1.2) 
[ C ifu = 0, 

where C 6 [liminf $(u),limsup 4>(it)]. 

Note that when j3{u) = u.\u\ m ~ l , m > 1, the partial differential equation (PDE) 
in (11.11) is nothing else but the classical porous media equation. In this case = 
| ii | 2 and in particular C = 0. 

Our main target is to analyze the case of an irregular coefficient ft. Indeed, we are 
particularly interested in the case when /3 is continuous excepted for a possible jump 
at one positive point, say u c > 0. A typical example is: 

P{u) = H(u - u c ).u, (1.3) 

H being the Heaviside function and u c will be called critical value or critical thresh- 
old. 

Definition 1.1. i) We will say that the PDE in dl.U . or /? is non-degenerate if there is 
a constant Co > such that $ > Co, on each compact of R + . 

ii) We will say that the PDE in dl.lt . or ft is degenerate if lim 4>(u) = 0. 

«->o+ 

Remark 1.2. i) We remark that f} is non-degenerate if and only if lim inf <&(u) > 0. 

u— S-0+ 

ii) We observe that j3 may be neither degenerate nor non-degenerate. 
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Of course, ft in dl -3b is degenerate. Equation (11.31) constitutes a model intervening in 
some self-organized criticality (often called SOC) phenomena, see |2) for a significant 
monography on the subject. We mention the interesting physical paper lfI31 . which 
makes reference to a system whose evolution is similar to the evolution of a "snow 
layer" under the influence of an "avalanche effect" which starts whenever the top of 
the layer is bigger than a critical value u c . 

We, in particular, refer to J9) (resp. 0), which concentrates on the avalanche phase 
and therefore investigates the problem dl.lt discussing existence, uniqueness and prob- 
abilistic representation when (3 is non-degenerate (resp. degenerate). The authors had 
in mind the singular PDE in dl.lt as a macroscopic model for which they gave a mi- 
croscopic view via a probabilistic representation provided by a non-linear stochastic 
differential equation (NLSDE); the stochastic equation is supposed to describe the 
evolution of a single point of the layer. The analytical assumptions formulated by the 
authors were Assumption(A) and the Assumption(B) below which postulates linear 
growth for /3. 

Assumption(B) 

There exists a constant c > such that < c\u\. 

Obviously we have, 

Assumption(B') 

There exists a constant c > such that |/?o(w)| < c\u\. Clearly (1 1.3b fulfills As- 
sumption (B). 

To the best of our knowledge the first author who considered a probabilistic repre- 
sentation (of the type studied in this paper) for the solutions of non linear deterministic 
partial differential equations was McKean 13T1 . However, in his case, the coefficients 
were smooth. From then on the literature steadily grew and nowadays there is a vast 
amount of contributions to the subject. 

A probabilistic interpretation of (11.11) when (3(u) = u.\u\ m ~ l , m > 1 was provided 
in J3). For the same /3, though the method could be adapted to the case where /3 is 
Lipschitz, in ||27l . the author studied the evolution problem dl.lt when the initial con- 
dition and the evolution takes values in the class of probability distribution functions 
on M. He studied both the probabilistic representation and the so-called propagation 
of chaos. 

At the level of probabilistic representation, under Assumptions(A) and (B), sup- 
posing that uq has a bounded density, J9) (resp. 10) proves existence and unique- 
ness (resp. existence) in law for (NLSDE). In the present work we are interested in 
some theoretical complements, but the main purpose consists in examining numerical 
implementations provided by (NLSDE), in comparison with numerical deterministic 
schemes appearing in one recent paper, see |[T71l . 

Let us now describe the principle of the probabilistic representation. The stochastic 
differential equation (in law) rendering the probabilistic representation is given by the 
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following (NLSDE): 



t 



Y + J ®(u(s,Y s ))dW s , 



< 



u(t,-) 
k w(O,0 





Law density of Y t , Vi > 0, 
uq Lawoflo, 



(1.4) 



where W is a classical Brownian motion. The solution of that equation may be visual- 
ized as a continuous process Y on some filtered probability space (£2, J 7 , (J-"t)t>o,P) 
equipped with an (J^^o-Brownian motion W. 

Until now, theoretical results about well-posedness (resp. existence) for dl.4t were 
established when (3 is non-degenerate (resp. possibly degenerate) and in the case when 
wo 6 (L l n^°°) Even if the present paper concentrates on numerical experi- 
ments, two theoretical contributions are performed when $ is continuous. 

• Initially our aim was to produce an algorithm which allows to start even with a 
measure or an unbounded function as intial condition. Unfortunately, up to now, our 
implementation techniques do not allow to treat this case. 

A first significant theoretical contribution is Theorem 12.91 which consists in fact 
in extending the probabilistic representation obtained by 10 to the case when uq £ 
L 1 (M), locally of bounded variation outside a discrete set of points. 

• A second contribution consists in showing in the non-degenerate case that the 
mollified version of PDE in dl.ll) is in fact equivalent to its probabilistic representation, 
even when the initial condition uq is a probability measure. This is done in Theorem 



The connection between (11.41) and dl.lt is then given by the following result. 

Proposition 1.3. Let us assume the existence of a solution Y for \1.4§ . Let u(t,-) be 
the law density ofYt, t > 0, that we suppose to exist. 

Then u : [0,T]xR-> M + provides a solution in the sense of distributions of Ul.li 
with uq = u(0, ■). 

The proof is well-known, but we recall here the basic argument for illustration pur- 
poses. 

Proof. Let (p £ C^°(R), Y be a solution of the problem (11.41) . We apply Ito's formula 
to (fi(Y) to obtain : 



<p(Y t ) = <p(Y ) + / V '(Y s )^(u(s,Y s ))dW s + - / p"(Y s )<t> 2 (u(s,Y s ))ds. 



E2 





Taking the expectation we get : 
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Using then integration by parts and the expression of ft, the expected result follows. 



In the literature there are several contributions about approximation of non-linear 
PDE's of parabolic type using a stochastic particles system, with study of the chaos 
propagation. We recall that the chaos propagation takes place if the components of a 
vector describing the interacting particle system become asymptotically independent, 
when the number of particles goes to infinity. Note that, physically motivated applica- 
tions can be found, for instance in numerical studies in hydro- or plasma-physics; |fl9l 
and 11231 are contributions expressing a heuristic or formal point of view. 

When the non-linearity is of the first order, a significant contribution was given by 
E3; iflOl [TT1 performed the rate of convergence, IT321 provided a chaos propagation 
result. We also quote lfl6l . where authors provided a propagation of chaos result for 
the Burger's equation. 

In the case of porous media type equation in (11.11) with ft Lipschitz, l28ll investigated 
the probabilistic representation for (11.11) and a mollified related equation. There, the 
authors provided a rigorous proof of propagation of chaos in the case of Lipschitz 
coefficients, see Proposition 2.3, Proposition 2.5 and Theorem 2.7 of 11281 . 

Outside the Lipschitz case, an alternative method for studying convergence was 
investigated by Il33l[34ll35l , whose limiting PDEs concerned a class of equations in- 
cluding the case ftiu) = u+u 2 , u > 0. In fact l35l computed the numerical solution 
of a viscous porous medium equation through a particle algorithm and studied the L 2 - 
convergence rate to the analytical solution. More recent papers concerning the chaos 
propagation when ft(u) = u 2 first and ft(u) = \u\ m ~ i u 1 m > 1 was proposed in 11381 
and EOll . 

As far as the coefficient ft is discontinuous, at our knowledge, up to now, there are 
no such results. As we announced, we are particularly interested in an empirical in- 
vestigation of the stochastic particle algorithm approaching the solution u of (11.11) at 
some instant t, in several situations with regular or irregular coefficient. We recall 
that u(t, •) is a probability density. That algorithm involves Euler schemes of stochas- 
tic differential equations, Monte-Carlo simulations expressing the empirical law and 
non-parametric density estimation of u(t, ■) using Gaussian kernels, see P51 for an 
introduction to the kernel method. This technique crucially depends on the window 
width e of the smoothing kernel. Classical statistical tools for choosing that parameter 
are described for instance in l45l . where the following formula for choosing the opti- 
mal bandwidth e, in the sense of minimizing the asymptotic mean integrated squared 
error (MISE), is given by 



where, n is the sample size and || • || denotes the classical Z 2 (M) norm. 

Of course, the above expression does not yield an immediately practicable method 
for choosing the optimal e since dl -5b depends on the second derivative of the density 



□ 




(1.5) 
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u, which we are trying indeed to estimate. Therefore, several techniques were pro- 
posed to get through this problem. First, a natural and easy approach, often called the 
rule of thumb, replaced the target density u at time t in the functional [jc^uH, by a 
reference distribution function. For instance, l45ll assumed that the unknown density 
is a standard normal function and obtained the following practically used formula 



a being the empirical standard deviation. A version which is more robust to outliers 
in the sample, consists in replacing a by a measure of spread of the variance involving 
the interquartile range. For instance, see l45ll for detailed computations. 

The oversmoothing methods rely on the fact that there is a simple upper bound for 
the MISE-optimal bandwidth. In fact, |48l . gave a lower bound for the functional 
|| it || and thus an upper bound for e in dl.51) : it proposed to use this upper bound as 
an optimal window width, see also ll49l for histograms. 

The two methods above seem to work well for unimodal densities. However, they 
lead to arbitrarily bad estimates of the bandwidth e, when for instance, the true density 
is far from being Gaussian, especially when it is a multimodal law. 

The least squares cross validation (LSCV) method aimed to estimate the bandwidth 
that minimizes the integrated squared error (ISE), based on a "leave-one-out" kernel 
density estimator, see l40l[T3l . The problem is that, for the same target distribution, 
the estimated bandwidth through different samples has a big variance, which produces 
instability. 

The biased cross-validation (BCV) approach, introduced in PD minimizes the 
score function obtained by replacing the functional HS^-itH in the formula of the MISE 
by an estimator ||S£ a .zt||, where u is the kernel estimator of u. In fact, RD proposed 
the use of the minimizer of that score function as optimal bandwidth. This method 
seemed to be more stable than the LSCV but still has large bias. The slow rate of 
convergence of both the LSCV and BCV approaches encouraged significant research 
on faster converging methods. 

A popular approach, commonly called plug-in method, makes use of an indirect 
estimator of the density functional ||c^ x u|| in formula ( 11.5b . This technique comes 
back to the early paper llBD : in this framework the estimator of j|<9^ x u|| requires the 
computation of a pilot bandwidth h, which is quite different from the window width 
£ used for the kernel density estimate. Indeed, this optimal bandwidth h depends on 
unknown density functionals involving partial derivatives greater than 2. Following an 
idea of BUI , one could express h iteratively through higher order derivatives. In this 
spirit, the natural associated problem consists in estimating for some positive integer 
s, the quantity ||9* s tt||, in terms of ||^^,u|| for some positive integer £; an i-stage 
direct plug-in approach may consist in replacing the norm ||9 s J^u|| by the norm of 
the s + £ derivative of a Gaussian density. In the present paper we implement this 




(1-6) 
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idea with s = i = 2. Important contributions to that topic were l42l and f26l who 
improved the method via the so-called "solve-the-equation" plug-in method. By this 
technique, the pilot bandwidth h used to estimate Hc^wH, is written as a function 
of the kernel bandwidth e. We shall describe in Section [4] in details this bandwidth 
selection procedure applied in the case of our probabilistic algorithm. 

We point out, that a more recent tool was developed in l52l which improved the idea 
in l42"l[26l in the sense that ll52l did not postulate any normal reference rule. However, 
the numerical experiments that we have performed using the Matlab routine developed 
by the first author of |52l have not produced better results in the case when is defined 
bv Oil. 

In the paper we examine empirically the stochastic particle algorithm for approach- 
ing the solution to the PDE in the case /3(u) = u 3 and in the case (3 given by dl -3b . For 
this more peculiar case, we compare the approximation with the one obtained by one 
recent analytic deterministic numerical method. 

Problems of the same type as dl.lt . in the case when (3 is Lipschitz but possibly 
degenerate, were extensively studied from both the theoretical and numerical deter- 
ministic points of view. In general, the numerical analysis of (11.11) is difficult for at 
least one reason: the appearance of singularities for compactly supported solutions in 
the case of an irregular initial condition. An usual technique to approximate (11.11) in- 
volves implicit discretization in time: it requires, at each time step, the discretization 
of a nonlinear elliptic problem. However, when dealing with nonlinear problems, one 
generally tries to linearize them in order to take advantage of efficient linear solvers. 
Linear approximation schemes based on the so-called non linear Chernoff 's formula 
with a suitable relaxation parameter and which arises in the theory of nonlinear semi 
groups, were studied for instance in JS). We also cite ||29l , where the authors approx- 
imated degenerate parabolic problems including those of porous media type. In fact, 
they used nonstandard semi-discretization in time and applied a Newton-like itera- 
tions to solve the corresponding elliptic problems. More recently, different approaches 
based on kinetic schemes for degenerate parabolic systems have been investigated in 
JO- Finally a new scheme based on the maximum principle and on the perturbation 
and regularization approach was proposed in 1*391 . 

At the best of our knowledge, up to now, there are no analytical methods dealing 
with the case when (3 is given by (11.3K However, we are interested in a sophisticated 
approach developed in IfTTl and which appears to be best suited to describe the evolu- 
tion of singularities and efficient for computing discontinuous solutions. In fact, IfTTl 
focuses onto diffusive relaxation schemes for the numerical approximation of nonlin- 
ear parabolic equations, see 12511241 . and references therein. Those relaxed schemes 
are based on a suitable semi-linear hyperbolic system with relaxation terms. Indeed, 
this reduction is carried out in order to obtain schemes that are easy to implement. 
Moreover, with this approach it is possible to improve such schemes by using differ- 
ent numerical approaches i.e. either finite volumes, finite differences or high order 
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accuracy methods. 

In particular, the authors in ffTTIl coupled ENO (Essentially Non Oscillatory) interpo- 
lating algorithms for space discretization, see 04], in order to deal with discontinuous 
solutions and prevent the onset of spurious oscillation, with IMEX (implicit explicit) 
Runge-Kutta schemes for time advancement, see l36ll . to obtain a high order method. 
We point out that ifTTll studied convergence and stability of the scheme only in the case 
when /3 is Lipschitz but possibly degenerate and uq £ L 1 (R). 

As a byproduct of numerical experiments we can forecast the longtime behavior of 
u(t, ■) where (t, x) h-> u(t, x) is the solution of the considered PDE. We can reasonably 
postulate that the closure of {u 6 L J (R), u > 0, J R u(x)dx = 1 | j3(u) = 0} is a 
limiting set, provided it is not empty as in the case /3(u) = u 3 . 

The paper is organized as follows. Section|2]is devoted to a survey of existence and 
uniqueness results for both the deterministic problem dl.lt and the non-linear SDE 
(ll.4t rendering the probabilistic representation of (ll.lt . We in particular, recall the 
results given by authors of J9j and we provide some additional results in the case 
when the initial condition of (ll.lt belongs to L l (R) but it is not necessarily bounded. 

In Section[3] we settle the theoretical basis for the implementation of our probabilis- 
tic algorithm. We first approximate the NLSDE (11.41) by a mollified version replacing 
u(t, ■), the law density of Y t , by a given smooth function. We then construct an in- 
teracting particle system for which we supposed that propagation of chaos result is 
verified. We drive the attention on Theorem 13 .21 which links the mollified PDE (13.31 ) 
with its probabilistic representation. 

Section |4] is devoted to the numerical procedure implementing the probabilistic al- 
gorithm. We first introduce an Euler scheme to obtain a discretized version of the 
interacting particles system defined in Section [3j We then discuss the optimal choice 
of the window width e. 

In Section [5] we describe the numerical deterministic approach we use to simulate 
solutions of (ll.lt . In fact, following ffTTIl . we first use finite differences and ENO 
schemes for the space discretization, then we perform an explicit Runge-Kutta scheme 
for time integration. 

In Section|6] we proceed to the validation of the algorithms. In fact, the first numer- 
ical experiments discussed in that section concern the classical porous media equation 
whose exact solution, in the case when the inial condition is a delta Dirac function, 
is explicitly given by the so-called Barrenblatt-Pattle density, see ff4|- Then, we con- 
centrate on the Heaviside case, i.e. with j3 of the form (ll.3t . In fact, we perform 
several test cases according to the critical threshold u c and to the initial condition uq. 
Finally, we conclude this section by some considerations about the longtime behavior 
of solutions of dl.lt . 
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2. Existence and uniqueness results 

We start with some basic analytical framework. If / : R — > R is a bounded function we 
will denote ||/||oo = sup \ f(x)\. By 5(R) we denote the space of rapidly decreasing 

infinitely differentiable functions cp : R ->■ R. We denote by M{R) and M+(R) the 
set of finite measures and positive finite measures respectively. 



2.1. The deterministic PDE 

Based on some clarifications of some classical papers |[6l[l4l[7]|, [9) states the following 
theorem about existence and uniqueness in the sense of distributions (in a proper way). 

Proposition 2.1. Let u e (l^fl-^ 00 ) (R), u > 0. We suppose the validity of 
Assumptions (A) and (B). Then there is a unique solution in the sense of distributions 



u 



e (L l f)L°°W,T}xm)of 



| d t u e \d 2 xx l3{u), 
y u(0,x) = u (x), 

in the sense that, there exists a unique couple (u,r} u ) 6 ((L 1 |~)L°°)([0, T] x R)) 2 
such that 

u(t,x)ip(x)dx = J uo(x)(p(x)dx + — J ds J r] u (s, x)ip"(x)dx, V</? 6 S'(R) 
and 

T]u(t, x) G (3(u(t, x)) for dt ® dx-a.e. (t, x) e [0, t] x R 

Furthermore, \\u(t, •) I loo < 1 1 w o| |oo /<"" every t € [0, T] an<i there is a unique version 
of u such that «£ C([0,T] jL^R)) (c ^([O.T] x R)). 

One significant difficulty of previous framework is that the coefficient ft is discon- 
tinuous; this forces us to consider /3 as a multivalued function even though u is single- 
valued. Being /3, in general, discontinuous it is difficult to imagine the level of space 
regularity of the solution u(t,-) at time t. In fact, Proposition 4.5 of Q says that al- 
most surely rj u (t, ■) belongs dt-a.e in ^(R) if u G (L 1 f| (R). This helps in 
some cases to visualize the behavior of u(t, ■). The proposition below makes some 
assertions when (3 is of the type of d 1 -3b . which constitutes our pattern situation. 

Proposition 2.2. Let us suppose uq 6 (L 1 p|L°°) (R) and (3 defined by (11.3b . For 

t > 0, we denote by 

E t — {x\ u(t,x) = u c }. 

For almost all t > 0, 
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( i) Et has a non empty interior; 

(ii) every point of E® is either a local minimum or a local maximum. 

Remark 2.3. The first point of the previous proposition means that at almost each time 
t > 0, the function u(t, ■) remains constant on some interval. 

The second point means that if the function u(t, ■) crosses the barrier u c , it has first 
to stay constant for some time. 

Proof of Proposition 12. 21 For the sake of simplicity we fix t > such that r) u (t, ■) £ 
H 1 (M.) and we write u = u(t, ■) , r} u = rj u (t, ■). 

(i) Since i] u £ H l (R) it is continuous, then the set Vq = {x £ R| r) u (x) £]0, u c [} 
is open. If rj u (x) £]0, u c [ necessarily we have u(x) = u c ; in fact, if u{x) < u c 
then r] u (x) = and if u(x) > u c then r) u (x) = u(x) > u c . Since Vq is open 
and it is included in E$ the result is established. 

(ii) Suppose the existence of sequences (x n ) and (y n ) such that with 
u(x n ) < u c and y n — >■ y with u(y n ) > u c . By continuity of r\ u we have 

Vu{x n ) =0^0 = t] u (x) 

n— >oo 

U(y n ) = Vu(Vn) ->■ Vu(x) = 0, 
n— >oo 

this is not possible because u(y n ) > u c for every n. 

□ 

If uq £ .M(R), we do not know any existence or uniqueness theorem for dl.ll) . Our 
first target consisted in providing some generalization to Proposition 12.11 in the case 
when uq is a finite measure. A solution in that case would be, n:]0,T]xl4L' (R) 
continuous and such that 

lim u(t, dx) = uo(dx), 
t->0 

weakly and where u(t, dx) denotes u(t, x)dx. This is still an object of further tech- 
nical investigations. For the moment, we are only able to consider the case uq having 
a L^R) density still denoted by uo, not necessarily bounded as in Proposition [2J] at 
least when $ characterized by dl.2t is continuous. In particular j3 is also continuous, 
but possibly degenerate. In that case, we can prove existence of a distributional solu- 
tion to dl.lt . Even though this is not a very deep observation, this will settle the basis of 
the corresponding probabilistic representation, completely unknown in the literature. 
In fact, we provide the following result. 

Proposition 2.4. Let uq £ L (R). Furthermore, we suppose that Assumption(A) and 
Assumption(B) are fulfilled. We assume that $ is continuous on R + . 
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(1) There is a solution u, in the sense of distributions, to the problem 

f d t u(t,x) = \dl x P{u{t,x)), t G [0,oo[, 
u(0,x) = uo(dx), 

in the sense that for every a G 5(M) 

u(t, x)a(x)dx = I uo(x)a(x)dx + — f ds f a" (x)f3(u(s,x))dx. 

(2.3) 

(2) If uq is locally of bounded variation excepted eventually on a discrete number 
of points Do, then <t>(u(t, ■)) has at most countable discontinuities for every 

te{0,T}. 

Proof (l)Let uo G L l (R), Uq = uo*4>±, N G N*, where 4> 

is a kernel with compact 

JV 

support and </>j_(x) = N(j)(Nx), x G R. So is of class C 1 , therefore locally 

with bounded variation. Since ||wg ||oo < \\4>± l|oo||wo||i,i then G (I/'H-k 
Moreover, we have 

/ \uq (x) — uo(x)\dx — >■ 0, as N — >■ +oo. 



On one hand, according to Proposition 12. 11 there is a unique solution u N of (12.31) . i.e. 

for every a G 



M JV (f,a;)a(a;)(ix = / v$ {x)a{x)dx + i / ds f a"(x)/3{u N (s,x))dx. (2.4) 
JR 2 Jo J K 

On the other hand, according to Corollary 8.2 in Chap IV of [43], we have 

sup / \u N (t, x) — u(t, x)\dx — > 0, as N — > +oo. (2.5) 

KTJl 

Therefore, there is a subsequence (Nk)keN such that 

u k (t, x) ^ u(t, x) dt ® drc-a.e., as fc — )■ +oo. 
Since /3 is continuous, it follows that 

(3(u Nk {t,x)) ->• 0{u(t,x)) dt^dx-a.e., as fe ->■ +oo. 
Consequently, ( 12.41 ) implies 

t(i, x)a(ar)da; = /" uo(a;)a(a;)da; + lim — f ds f a" (x)fl(u Nk (s,x))dx. 

JR k^+oo 2 Jo J R 

(2.6) 
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In order to show that u solves (12.31 ) . we verify 



N( n _ / a„ I 



lim / ds / a"{x)P(u N {s,x))dx = / ds / a"(x)P(u(s, x))dx, (2.7) 
where for notational simplicity we have replaced with N. So, we can suppose that 



N , „, o(„.N 



u, P{u n ) /3(u), eft ® dx-a.e. as TV -> +oo. (2.8) 



Since 1/3(^)1 < c\u N \ mdu N -> uinL^^T] x M), it follows that f3{u N ) are equi- 
integrable. Consequently, by (f2~8T ). -> /3(u) in ^'([O,^ x M), and therefore 

(12.71) follows. Finally, u solves equation ( 12.3I I. 

(2) For this purpose we state a lemma concerning an elliptic equation whose first 
statement item constitutes the kernel of the proof of Proposition 12. II 

Given / : M ->■ R, for /i 6 E, we denote 

/ /l (a;) = /( a; + / l )-/(x). 

Lemma 2.5. Let f e L 1 , X > 0. 

(i) There is a unique solution in the sense of distributions of 

u - AG9(u))" = /. 
f »J Lef X^ e a smooth function with compact support. Then for each h 

X (x)\u h (x)\dx < [ x (a)|/ h (x)|do: + CA|/i|||u|| L i, (2.9) 



where C is a constant depending on /3 and \- 

Proof of Lemma |231 (i) is stated in Theorem 4.1 of ||6l and Theorem 1 of Q. 

(ii) The statement appears in Lemma 3.6 of J3) in the case when / 6 L l f] L°° but 
the proof remains the same for / 6 L . □ 

We go on with the proof of Proposition l2.4l point (2). Let x be a smooth nonnegative 
function with compact support on M\Dq. We prove in fact 

limsup-^- J x( x )\ uh {ti x)\dx < ||uoAl|var + C f \u(s,x)\dsdx, (2.10) 

where || ■ || va r denotes the total variation and C is a generic universal constant. For 
this purpose, we proceed exactly as in the proof of Proposition 4.20 of Q making 
use of Lemma 12.51 Inequality ( 12.101 1 allows, similarly as in J3) to show that u(t, ■) 
restricted to any compact interval of M\Dq has bounded variation. Therefore it has at 
most countable discontinuities. Consequently, &(u(t, ■)) has the same property since 
is supposed to be continuous. □ 
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2.2. The non-linear stochastic differential equation (NLSDE) 

Definition 2.6. We say that a process Y is a solution to the NLSDE associated to 
problem dl.ll . if there exists \ belonging to L°°([0, T] x R) such that; 

Yt = Y + tix(s,Y s )dW s , 
X(t,x) G &(u(t,x)), for dt ® dx - a.e. (t,x) G [0, T] x R, (2 11) 
u(t,x) = Law density of Y t , t > 0, 
u(0, •) = uo, 

where W is a Brownian motion on some suitable filtered probability space (£2, J 7 , ( .Ft )t>o 
In particular, the first identity of (12.111 1 holds in law. We introduce a notion appearing 
in El- 
Definition 2.7. We say that ft is strictly increasing after some zero if there is a constant 
c > such that 

i) /% c] =0. 

ii) /3 is strictly increasing on [c, +oo [. 

Up to now, two results are available concerning existence and uniqueness of solu- 
tions to (12.111) . In fact, the first one is stated in the case where ft is not degenerate and 
the second one in the case when ft is degenerate, see respectively J9][3l. We summarize 
these two results in the following theorem for easy reference later on. 

Theorem 2.8. Let uq g L l f*\L°° such that uq > Oand J K uo(x)dx = 1. Furthermore, 
we suppose that Assumptions (A) and (B) are fulfilled. 

(i) If ft is non-degenerate then it exists a solution Y to d2.1U . unique in law. 

( ii) Suppose ft is degenerate and either ft is strictly increasing after some zero or uo 
has locally bounded variation. Then there is a solution Y not necessarily unique 
to (PTTTT) . 

A step forward is constituted by the proposition below. This provides an existence 
result for the NLSDE, when uq is not necessarily bounded at least whenever <& is 
continuous. This does not require a non-degenerate hypothesis on ft. 

Theorem 2.9. Let uq G L'(R) having locally bounded variation except on a discrete 
set of points Dq. Furthermore we suppose that Assumption(A) and Assumption(B) are 
fulfilled. We assume that <E> is continuous on R + . 

The probabilistic representation related to dl.ll) holds, i.e. there is a process Y 
solving (1 1 -4b in law. 
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Proof. Let Uq be the function considered at the beginning of the proof of Proposition 
2.41 According to Theorem l2.8l let Yq be the solution to 

y t N = Y Q N + J*<t>(u N (s,Y s N ))dW s , 
u N (t,-) = Law density of Y t N , V t > 0, (2.12) 
u N (0,-) = u». 

Since <I> is bounded, using Burkholder-Davies-Gundy inequality one obtains 

E(Y t N -Yf) 4 < const(t - s) 2 . 

This implies ( see for instance Problem 4.11, Section 2.4 of f30l ) that the laws of 
Y N ,N > 1 are tight. Consequently, there is a subsequence Y k := Y Nk converging 
in law (as C([0, T])-valued random elements) to some process Y. We set u k = u Nk , 

where we recall that u k (t, ■) is the law of Y t k . We also set X k = Y k — Y^. Since 

t 

[X k ] t = J $> 2 (u k (s, Y k ))ds and <1> is bounded, the continuous local martingales X k 

o 

are indeed martingales. 

By Skorohod's theorem there is a new probability space (£2, J ', P) and processes 
Y k , with the same distribution as Y k so that Y k converges P-a.s. to some process 
Y, of course distributed as Y, as C([0, T])-valued random element. In particular, the 
processes X k = Y k — Y^ remain martingales with respect to the filtration generated 
by Y k . We denote the sequence Y k (resp. Y), again by Y k (resp. Y). 

We now aim to prove that 

Y t = Y + f ®(u(s,Y s ))dW s , (2.13) 

for some standard Brownian motion W with respect with some filtration (J 7 *). 

We consider the stochastic process X (vanishing at zero) defined by Xt = Yt — Y$. 
We also set again X k = Y k — Yq. Taking into account Theorem 4.2 in Chap 3 of l30l . 
to establish ( I2.13I I. it will be enough to prove that X is an ^-martingale with quadratic 
variation [X] t = J* * <& 2 (u(s, Y s ))ds, where y is the canonical filtration associated with 
Y. 

Let s, t G [0, T], with t > s and ip a bounded continuous function from C([0, s]) to 
R. In order to prove the martingale property for X, we need to show that 

E[(X t -X a )1>(Y r ,r<s)] = 0. (2.14) 

Since Y k are martingales, we have 



E {X k -X k )^{Y k ,r<s) 



0. (2.15) 
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Consequently ( 12.141 1 follows from ( 12.151) and the fact that Y k -t Y a.s. (X k X 
a.s.) as C([0, T])-valued random process. In fact for each t > 0, X t fc -> X t in L 1 (H) 
since (X t fc , fc e N) is bounded in L 2 (£l). 

It remains to show that X 2 — J^^> 2 (u(s,Y s ))ds, t £ [0, T], defines an -martingale, 
that is, we need to verify 



E 



Xt ~ X 2 S - <5 2 (>(r, Y r ))drj ^{Y r ,r < 



0. 



We proceed similarly as in the proof of Theorem 4.3 in J9) but even with some simpli- 
fication. For the comfort of the reader we give a complete proof. 
The left-hand side decomposes into I 1 (k) + 1 (k) + 1 (k), where 



I\k) 



I\k) 
I\k) 



E 
E 

E 

E 



X t — X e 



® 2 (u(r,Y r ))drj ip(Y r ,r < s) 
(X*) 2 - (X k ) 2 - j\ 2 {u{r,Y k ))dr^ ^{Y k ,r < s) 



{X k f-{X k ) 2 - I &{u"{rX))dr)i>{Y r «,r< s ) 



{<P 2 {u k (r, Y k )) - <D z (w(r, Yf)))dr ) ^(Y*, r < s) 



We start by showing the convergence of I 3 (k). Now, tp(Y k , r < s) is dominated by a 
constant C. Clearly we have 

l\k)<C f dr [ \<t> 2 (u k (r,y))-<i> 2 (u(r,y))\u k (r,y)dy. 

The right hand side of this inequality is equal to C[J l (k) + J {k)\, where 

J l (k) = J dr Jj<P 2 (u k (r,y))-<P 2 (u(r,y))\(u k {r,y)-u(r,y)^dy, 

J\k) = [ dr [ \<Z> 2 (u k (r,y))-<t> 2 (u(r,y))\u(r,y)dy. 



Since u k -> u in C([0,T];L l ) and <I> 2 is bounded then lim J l {k) = 0. 
Furthermore, there is a subsequence (fc„) ne pj such that 



u kn (r, y) — > u(r, y) dr ® dy — a.e. as n-> +oo. 
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Since <I> 2 is continuous, it follows that 

(S> 2 (u kn (r,y)\ ->■ <I> 2 (u(r,y)) dr ® dy - a.e., TV ->• +00 
On the other hand, since 

|$ 2 (u k "(r,y)) -<t> 2 (u(r,y)) \ < 2 sup <t> 2 (u)\u(r, y)\, 

Lebesgue's dominated convergence Theorem implies that lim J 2 (k) = 0. 

k— >+oo 

Now we go on with the analysis of I 2 (k) and I x (k). I 2 (k) equals zero since X k is 

i 

a martingale with quadratic variation given by [X] t = J <& 2 (u k (r, Y k ))dr. 



Finally, we treat I x (k). We recall that X k — > X a.s. as a random element in 
C([0,T}) and that the sequence E((X k ) 4 ) is bounded, so (X k ) 2 are uniformly inte- 
grable. 

Therefore, we have 



E \{{X t ) 2 - (X s ) 2 ) Tp(Y r ,r < s)] - E [([X k f - {X k f) ^(Y k ,r < s) 
when k goes to infinity. It remains to prove that 

(V(u(r, Y r )) - <t> 2 (u(r, Y k ))^j $ {Y r ,r < s) dr 







E 



->0. (2.16) 



Now, for fixed <ir-a.e., r G [0, T], the set S(r) of discontinuities of <J>(-u(r, .)) is 

countable because of Proposition 12.41 point (2). The law of Y r has a density and it is 

therefore non-atomic. Let N(r) be the event of all to £ D, such that Y r {uS) belongs to 

S{r). The probability of N{r) equals E(l smO^)) = / K l s(r){y)dv(y) = 0, where v 

is the law of Y r . Consequently N(r) is a negligible set. 

For uj 4 N (r), we have lim <I> 2 (u(r, Y k (uj))) = <S> 2 (u(r, Y r (ui))). Since ^> is 
k^+00 

bounded, Lebesgue's dominated convergence theorem implies (12.161) . 

Concerning the question wether u(t, .) is the law of Yt, we recall that for all t, Y k 
converges (even in probability) to Y t and u k (t, .), which is the law density of Y k , goes 
to u(t, .) in L J (]R). By the uniqueness of the limit in (12.31) , this obviously implies that 
u(t, .) is the law density of Y t . □ 



3. Some complements related to the NLSDE 
3.1. A mollified version 

We suppose here that uo(dx) is a probability measure. Let Yq be a random variable 
distributed according to u (dx) and independent of the Brownian motion W. 
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In preparation to numerical probability simulations, we define K £ for every e > 0, 
as a smooth regularization kernel obtained from a fixed probability density function K 
by the scaling : 

KJx) = -K (-) , x e R. (3.1) 

We suppose in this section that <I> is single valued, therefore continuous. This hypoth- 
esis will not be in force in Sections 0] and [6] 

In this subsection we wish to comment about the mollified version of the NLSDE 
(11.41 ). given by 

Y{ = Yo + f®((K £ *ve)( S ,Y s e))dW s , 
* v £ {t,-) = Law of Yf, Vt > 0, (3 ' 2) 

and its relation to the nonlinear integro-differential PDE 

d t v £ {t,x) = \d 2 xx (® 2 (K £ *v £ (t,x))v e (t,x)) , (t,x) e ]0,+oo[ x R, 
v £ (0,-) = u . 

(3.3) 

where, t h-> u E (t, ■) may be measure-valued. 

Remark 3.1. (i) When $ is Lipschitz, the authors of l28l proved in Proposition 
2.2, that the problem d3.2l l is well-posed. Their proof is based on a fixed point 
theorem with respect to the Kantorovitch-Rubeinstein metric. 

(ii) At our knowledge, there are no existence and uniqueness results for (13.31) at least 
when <I> is not smooth. 

(iii) By Ito's formula, similarly to the proof of Proposition 1 1.31 it is easy to see that a 
solution Y £ of (13.21) provides a solution v £ of ( 13.31 . in the sense of distributions. 

When (3 is non-degenerate it is possible to show that formulations d3.3l l and d3.2l l 
are equivalent. In particular we have the following result. 

Theorem 3.2. We suppose that (3 is non-degenerate and e > is fixed. 

(1) IfY £ isa solution of i3.2\ then v e : [0, T] -> M(R), where v £ (t,-) is the law 
ofY t £ , is a solution of (13.31) and fulfills the following property 

(P) v £ has a density, still denoted v £ such that: (t, x) i— > v £ (t, x) G Z/ 2 ( [0, T] xR). 

(2) If v £ is a solution to d3.3l) fulfilling (P) then there is a process Y = Y £ solving 
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Proof. (1) If Y E is a solution to ClD by Remark|3jJ(iii) it follows that v e fulfills (f33T >. 
On the other hand, since K 6 * v £ is bounded and $ is lower bounded by a constant 

C e on [— inf K E * w e , sup K £ * v e ] it follows that a(t, x) = 4> 2 {K £ * x)) is lower 
bounded by C E . 

Using then Exercise 7.3.3 of ll46l . i.e., Krylov estimates, it follows that for every 
smooth function / : [0, T] xl^l with compact support, we have 



eU f(Y s £ )ds\ <const||/|| iW]; 



Then, developing the left hand side to obtain 

rT 

l,T]xI 



/ ds [ f(y)v £ (s,y)dy <const\\f\\ L in Q 7 
Jo Jr 



we deduce that (P) is verified. 

(2) We retrieve here some arguments used in the proof of Proposition 4.2 of Q. 
Given v = v £ , by Remark 4.3 of (9), see also Exercise 7.3.2-7.3.4 of 1146*1 . we can 
construct a unique solution Y = Y e in law to the SDE constituted by 

i 

Y t = Y + J a{ Sl Y s )dW s , (3.4) 
o 

where here a(t, x) = 4> 2 (K e * v(t, x)) . Indeed, this is possible again because a is 
Borel bounded and lower bounded by a strictly postive constant. 
A further use of Ito's formula says that the law z(t, dx) of Y t solves 




\d 2 xx (a(t, .)z(t, .)) 
uo, 



(3.5) 



in the sense of distributions. 

Using again Krylov estimates as in the second part of the proof of point (1), it 
follows that z admits a density (t, y) n> pt(y) which verifies p 6 Z 2 ([0, T] x M). 
This shows that Hypothesis 13 .41 in Theorem 13.31 below is fulfilled, which implies that 

v = z. n 

Theorem [33] was stated and proved in |9|, see Theorem 3.8. 

Theorem 3.3. Let a be a Borel nonnegative bounded function on [0, T] x M. 

Let Zi : [0, T] — >■ JA + (M.), i = 1,2, be continuous with respect to the weak topology 
on finite measures on Ai(R). 
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Let z° be an element of M + (WL). Suppose that both z\ and zi solve the problem 
dtz = d\ x [az) in the sense of distributions with initial condition z(0, ■) = z . 
More precisely 



for every t G [0, T] and any 6 C£°(K). 

Then (zi — Z2)(t,-) is identically zero for every t, if z := z\ — Z% satisfies the 



Hypothesis 3.4. There is p : [0, T] x R -> R belonging to L 2 ([k, T] x M) for every 
k > such that pit, ■) is the density of z(t, ■) for almost all t g]0, T]. 

3.2. The interacting particles system 

We recall that in this paper, we want to approximate solutions of problem dl.U . For 
this purpose we will concentrate on a probabilistic particles system of the same nature 
as in f28ll when the coefficients are Lipschitz. 

In general, the particles probabilistic algorithms for non linear PDEs are based on 
the simulation of particles trajectories animated by a random motion. The solution 
of the PDE is approximated through the smoothing of the empirical measure of the 
particles, which is a linear combination of Dirac masses centered on particles positions. 
This procedure is heuristically justified by the chaos propagation phenomenon which 
will be explained in the sequel. 

The dynamics of the particles is described by the following stochastic differential 
system: 



where W = (W , . . . , W n ) is an n-dimensional Brownian motion, (^o l )i<i<n is a 
sequence of independent random variables with law density uq and independent of the 
Brownian motion W and K £ is the same kernel as in Subsection[3jl. 

Remark 3.5. If <I> in the system of ordinary SDEs (13.61 1 were not continuous but only 
measurable, that problem would not have necessarily a solution, even if /3 were non- 
degenerate. In fact, contrarily to (13.41) . here n > 2. Since <I> is continuous, then (13.61 1 
has at least a solution; if $ is non-degenerate even uniqueness holds, see Chapter 6 
and 7 of Egl . 

Now, owing to the interacting kernel K £ , the particles motions are a priori depen- 
dent. For a given integer n, we consider (Y^ ,e,n , . . . , Y™' 6 ' 71 ) as the solution of the 




following: 




(3.6) 
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interacting particle system ( 13.61) . Propagation of chaos for the mollified equation hap- 
pens if for any integer m, the vector (Y^ 1,e '", . . . ,l^ m ' e ' n )„> m converges in law to 
fj>t® m where /^ is the law of Yf the solution of d3.2l l. 

A consequence of chaos propagation is that one expects that the empirical measure 

of the particles, i.e. the linear combination of Dirac masses denoted = 7t S $Y j ' £ ' n 

i=i ' 

converges in law as a random measure to the deterministic solution v E (t, .) of the reg- 
ularized PDE (13.31) which in fact depends on e. This fact was established for instance 
when /3 is Lipschitz, in Proposition 2.2 of fl28ll . On the other hand when e goes to 
zero, the same authors show that v £ converge to the solution u of dl.ll) . They prove 
the existence of a sequence (e(n)) slowly converging to zero when n goes to infinity 

1 n 

such that the empirical measure ^ Yl ^x-j,". converges in law to u, see Theorem 

j=i t 

2.7 of l28ll . One consequence of the slow convergence is that the regularized empirical 
measure 

1 ™ 

-l^ K e(n)(--Y t ) 

i=i 

also converges to u. Consequently, this probabilistic interpretation provides an algo- 
rithm allowing to solve numerically dl.ll ). 

We recall however that one of the significant object of this paper is the numerical 
implementation related to the case when, /3 is possibly discontinuous; for the moment 
we do not have convergence results but we implement the same type of algorithm and 
we compare with some existing deterministic schemes. 



4. About probabilistic numerical implementations 

In this section we will try to construct an approximation method for solutions u of 
(11.11) . based upon the time discretization of the system (13.61) . For now on, the number 
n of particles is fixed. 

In fact, to get a simulation procedure for a trajectory of each (Y?' £ ' n ), i = 1, . . . , n, 
we discretize in time: for fixed T > 0, we choose a time step At > and TV e N, such 
that T = NAt. We denote by = kAt, the discretization times for k = 0, . . . , N. 

The Euler explicit scheme of order one, leads then to the following discrete time 
system, i.e., for every i = 1, . . . ,n 

where at each time step t^, we approximate u(tk, .) by the smoothed empirical measure 
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of the particles : 

1 - 

l {t k ,x) = -Y j K £ {x-Xl k ), k = l,...,N, xeR, (4.2) 

71 * K 



u e,n ( 



n 



at each time step and for every i — 1, . . . , n, the Brownian increment [Wl k +1 ~~ ^t k 
is given by the simulation of the realization of a Gaussian random variable of law 
M(0,At). 

One difficult issue concerns the smoothing parameter e related to the kernel K e . It 
will be chosen according to the kernel density estimation. 

In fact from now on we will assume that K, as defined in (13.11) . is a Gaussian 
probability density function with mean and unit standard deviation. In this case, in 
(14.21) . the function u £,n {tk, ■) becomes the so-called Gaussian kernel density estimator 
of u(fk, ■) for every time step with k = 1, . . . ,N. 

Finally, the only unknown parameter in ( 14.21) . is e; most of the authors refer to it as 
the bandwidth or the window width. 

The optimal choice of e was the object of an enormous amount of research, because 
its value strongly determines the performance of u £ ' n as an estimator of u depends, 
see, e.g. l45l and references therein. The most widely used criterion of performance 
for the estimator ( 14.21 ) is the Mean Integrated Squared Error (MISE), defined by 

W$E{u e > n (t,x)} = E u J[u^ n (t iy )-u(t,y)} 2 dy 

( 

E„ [u £ > n (t, y)] - u(t, y)\ dy+ I V„ [u e '*(i, y)) dy, 



\ point-wise bias 



integrated point-wise variance 



where, E M and Y u are respectively the expectation and the variance of X\ , j = 1 , . . , n 
under the assumption that they are independent and distributed as u(t, ■). 

We emphasize that the MISE expression is the sum of two components: the inte- 
grated bias and variance. 

The asymptotic properties of ( 14.21 ) under the MISE criterion are well-known (see ll45l , ll52l "). 
but we summarize them below for convenience of the reader. 

Theorem 4.1. (Properties of the Gaussian kernel estimator) 

Under the assumption that e depends on n such that lim e = 0, lim ne = +oo 

n— »+oo n— >+oc 

and d\ x u is a continuous square integrable function, the estimator ( 14.21 ) has integrated 
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squared bias and integrated variance given by 

\\E u [u £ ' n (t,.) -u{t,.)] || 2 = -e 4 \\8l x u\\ 2 + o(e 2 ), n -> +00, (4.3) 

[v u [u^ n (t,y)]dy = —!— + o((ne)- [ ),n^+^. (4.4) 
J 2enV7r 

Remark 4.2. (i) Here | . || denotes the standard L 2 norm. The first order asymptotic 
approximation of MISE, denoted AMISE, is thus given by 

AMlSE{u E ' n (t,x)} = -e 4 \\d 2 x u(t,x)f + {2en^)~ l . (4.5) 

(ii) The asymptotically optimal value of e is the minimizer of AMISE and by simple 
calculus it can be shown (see j37ll . Lemma 4A) to be equal to e° pt defined in 
formula dl.5t . 

As argued in the introduction, we have chosen to use the "solve-the-equation" band- 
width selection plug-in procedure developed in 14211261 . to perform the optimal win- 
dow width of the Gaussian kernel density estimator u £,n of u, defined in d4.2l l. 



Remark 4.3. According to 11421 . for every positive integer s, the identity 

||aX^.)|| 2 = (-l) S / d 2 J sU (t,x).u(t,x)dx, 

Jr 

suggests the following estimator for that density functional: 

n £ i=\ j=\ v 6 / 

where, is the r th derivative of the Gaussian kernel K and d x rU is the r th partial 
spacial derivative of u. 

Inspired, by (II -5b . the authors of Il42ll26ll look for an approached optimal bandwidth 
for the AMISE as the solution of the equation 



(2n^||C« 7(E() '"(^)l| 2 )" 1/5 (4-7) 



where, ||^ x w 7 "'"'|| 2 is an estimate of ||^ x u|| 2 using (14.61) . for s = 2, and the pilot 
bandwidth 7(e), which depends on the kernel bandwidth e. The pilot bandwidth 7(e) 
is then chosen through an intermediate step which consists in obtaining a quantity h t 
minimizing the asymptotic mean squared error (AMSE) for the estimation of He^uH 2 . 
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AMSE is in fact some approximation via Taylor expansion of the MSE which is de- 
fined as follows: 

2 



MSE{||^y*'1 2 }=E 



2 „Ami2 



\d xx u 



\\d xx u\ 



Similarly one can define an analogous quantity for the third derivative 

2 



MSE{||d 3 3 u M || 2 } =E 



,3 n ,h,n\\2 



\d\u 



(4.8) 



(4.9) 



and related AMSE. Exhaustive details concerning those computations are given in 
BUI . In fact, the authors in ll50l computed those minimizers and provided the follow- 
ing explicit formulae 



hi 



n\\dl 3 u(t,x)\ 



1/7 



ht 



-2K®(0) 
n\\d 4 x4 u(t,x)\\ : 



1/9 



(4.10) 



where ht and h^ minimize the AMSE corresponding respectively to (14.81) and i 

Solving dl.51 . with respect to n and replacing n in the first equality of d4.10l l. gives 
the following expression of ht in term of et 



h t 

This suggests to define 

list) 



1 V7 



\\9 3 x Mt,x)\\ 2 
4^K< 4 \0)\\d 2 xx u h t> n (t,x)\\ 2 



-5/7 



1/7 



-5/7 



(4.11) 



where, ||<9 2 , E w' l t' n || 2 and ||5 3 ,u 
bandy 

2iv~ (4) (0) 



t ,nn2 me estimators of Hc^till 2 and |<9 3 3 -u|| 2 using 



formula (14.61 1 and pilot bandwidths h\ and h\ given by 

1/7 



-2iv~( 6) (0) 
n\\dtMt7x)f 



1/9 



||d 3 3 ii(i, x)|| 2 and || d 4 4 ii(i, a;) || 2 will be suitably defined below. Indeed, h\ and /i 2 
estimate h t and denned in (14. 10b . 

According to the strategy in l42l l50l . we will first suppose that d \u(t, x) and 
d^u(t, x) are the third and fourth partial space derivatives of a Gaussian density with 
standard deviation a t of X t . In a second step we replace at with the empirical standard 
deviation &t of the sample X\, . . . , X™. This leads naturally to 
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Coming back to (14.71) . where 7(st) is defined through (14.111) . it suffices then to perform 
a root-finding algorithm for it at each discrete time step tk, in order to obtain the 
approached optimal bandwidth s tk . 



5. Deterministic numerical approach 

We recall that the final aim of our work is to approximate solutions of a nonlinear 
problem given by 

d t u(t,x) e ±dl x p(u(t,x)), t e [0,+oo[, 
u(0, x) = uq(x), 

in the case where ft is given by (11.3I ). Despite the fact that, up to now at our knowledge, 
there are no analytical approaches dealing such issues, we got interested into a recent 
method, proposed in IfTTl . Actually, we are heavily inspired by ifTTl to implement a 
deterministic procedure simulating solutions of (15.11 1 which will be compared to the 
probabilistic one. IfTTl handles with the propagation of a discontinuous solution, even 
though coefficient ft is Lipschitz. It seems to us that in the numerical analysis literature, 
IfTTl is the closest one to our spirit. We describe now the fully discrete scheme we will 
use for this purpose. 



5.1. Relaxation approximation 



The schemes proposed in 11171 follow the same idea as the well4inown relaxation 
schemes for hyperbolic conservation laws, see lf25l for a review of the subject. For 
the convenience of the reader, we retrieve here some arguments of IfTTl . where we 
recall that the coefficient ft is Lipschitz. In that case G, of course, becomes =. 

The equation (15.11 ) can be formally expressed by the first order system on R + x R : 




(5.2) 

v + \d x ft(u) = 0. 

(15.21) . is relaxed with the help of a parameter e > 0, in order to obtain the following 
scheme 

(d t u + d x v = 0, 

\d t v + ^d x ft(u) = -\v. 

Then, another function w : M + x R — >■ R is introduced in order to remove the non- 
linear term in the second line of system (15.31) . So, we obtain 



(5.4) 

ft(u)). 
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Note that d5.4l l is a particular case of the BGK system previously studied in lfl2l . In 
fact, authors of lfl2l proved that w (resp. v) converges to (3(u) (resp. — ^d x (3(u)), as 
e — > + . Furthermore, they showed the convergence of solutions of (15.41) to those of 
PDE (15.11) , in L l (M), as e goes to zero. 

Finally, we introduce a supplementary parameter if > 0, according to usual nu- 
merical analysis techniques; while preserving the hyperbolic character of the system. 
Therefore, we get 



d t u + d x v = 0, 

d t v + <f 2 d x w = -\v + (ip 2 - j£ 
d t w + d x v = -\{w - f3(u)). 



)d x w, 



(5.5) 



Now, setting 



u \ 

V 

w 



Hz) 



1 \ 

if 2 
10/ 



and g(z) 



the system (15.51 ) is rewritten in matrix form as follows 





-v + (if 2 £ - \)d x w 
/3(u) — w 



d t z + d x T{z) = -g(z). 

E 



Using the change of variable Z = P L z, where 



2ip 

\l -I J 



o =! A 

up 2 



and 



<p 





\ 





-<p 











o ) 



(5.6) 



we obtain 




• , , V + fW -v + <pw 
with U = — 7 r J — , V = „ 

2yj 2cp 



, W = u — w, 



(5.7) 



where, U, V, W are called characteristic variables. 
Since z = FZ, equation (15.61 ) leads to 



d t z + : 



l g(¥Z). 



(5.8) 
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By rewriting the system d5.8l l in terms of the characteristic variables, we obtain 

/ d t U + <pdJA \ 



d t V - (pd x V 



\ 



l g(¥Z). 



(5.9) 



Finally, solving d5.5l l is equivalent to the resolution of a three advection equations 
system, ( 15.91 1, with respectively a positive, a negative and a zero advection velocity. 

Remark 5.1. Note that we can deduce from ( 15.71) the following relation 

U = U + V + W. (5.10) 

5.2. Space discretization 

In the sequel of this chapter and in Annex [7] given two integers i < j, [[i, j]] , will 
denote the integer interval {i, i + 1, . . . , j}. We will now provide a space discretization 
scheme for system ( 15.91) . Let us introduce a uniform grid on [a, b] C R. 

We denote Xi = a — 4p + iAx , i e [[l,iV ffi ]] and x i+l / 2 = a + iAa; , i e [[0, N x ]} , 
where Ax = - a is the grid spacing and N x the number of cells. Note that Xj is the 

center of the interval [xi_i/ 2 , x i+ i/ 2 ]- Moreover, we denote the boundary conditions 
by u(t, a) = u a (t) and u(t, b) = u\,{t), for every t > 0. 

Then, we evaluate ( 15.91 ) on the grid of discrete points (xi) getting, 



Gi^ari), V* >0, Vie [[1,7V,]], 
(? 2 (*,a?i), Vi >0, Vie [[1,7V,]], 
G 3 (t,ari), V* >0, Vi e [[1, TV,]], 



(5.11) 



where, 



(Gi, G2, G-sY 



l g(¥Z). 



(5.12) 



Remark 5.2. We can easily deduce from ( 15.121 ). that for every t e]0, +00 [ and every 

3 

i e [1, JVal, we have : £ Gj(i, a?i) = 0. 



In order to ensure the convergence of the semi-discrete scheme (15.111 ) it is necessary 
to write it in a conservative form. To this aim, following ifTTl . we suppose the existence 
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of functions U and V such that 



U ^ x ) = ^ I U{t,y)dy, Vxe]a,b[, Vi > 0, 

x—Ax/2 
x-\-Ax/2 

V(t, x) = ^ / V(i, e]a, &[, V* > 0. 

x-Az/2 



n:+A:r/2 



Substituting in (15.111) . we obtain for every t > and every i 6 [1, iV x J, 
x<) + ^ x i+l/2 ) - U(t, Xi _ l/2 )) = Gi (*, xO, 



hx 



V(t,x i+l/2 ) - V{t,Xi_ 1/2 ) 



G 2 (t,xi), 



G3(t,Xi). 



(5.13) 



Let us now denote by Ui + \/ 2 {t) and Vj+i/2(i) the so-called semi-discrete numerical 
fluxes that approximate respectively U(t,Xi + \/ 2 ) and V(t, Xj+irt). For the sake of 
simplicity, we chose to expose only the calculations necessary to obtain the first semi- 
discrete flux U i+ \/ 2 (t), the same procedure being applied for the other one. 

In order to compute the numerical flux U i+ \/ 2 {t), we reconstruct boundary extrap- 
olated data U^^ 2 (t), from the point values Ui(t) = U(t,Xi) of the variables at 
the center of the cells, with an essentially non oscillatory interpolation (ENO) method. 
The ENO technique allows to better localize discontinuities and fronts that may appear 
when /3 is possibly degenerate; see Il22ll44ll for an extensive presentation of the subject. 
In fact, Wy +1 j 2 {t) (resp. Wr^jy 2 (t)) is calculated from an interpolating polynomial of 
degree d, on the interval [2^+1/2,2^+3/2] (resp. [2^-1/2,2^+1/2]) using a so-called ENO 
stencil, see ll44l and formula (17.5b in Annex 1770 

Next, we shall apply a numerical flux to these boundary extrapolated data. In order 
to minimize the numerical viscosity and according to authors of ffTTl . we choose the 
so-called Godunov flux, $q, associated to the advection equation 



and defined as follows 



d t u + d x f{u) = 0, 



min /(£), if a < 7, 

Q<4<7 



£ G [a, 7] = < 

where /(£) = ip£, with cp > 0. So we have, JJg [ a > 7] = V" 



max /(£), if 7 < a. 

. l<t<u 
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In fact, we set 

Vt>0, U i+y2 (t)=$ G [U7 +l/2 (t),U+ +l/2 (t)}. (5.14) 
Therefore, we obtain the following semi-discrete flux 

Vt>0, U i+l/2 (t) = ipU7 +1/2 (t). (5.15) 

Applying the previous procedure to compute V i+ i/ 2 (0 and replacing in (15.131) . we get 
for every t > and every i 6 [[1 , -/V x ]] , 

< f ^ Xi )-^{vt +1/2 {t)-Vt l/2 {t)) = G 2 (t, Xi ), (5.16) 

^(t,ar<) = G 3 (t, Xi ). 

Consequently summing up the three equation lines in (15.161) and using Remarks 15. II 
and 15.21 we obtain 

f t (t, Xi ) + £ (ur +l/2 (t)-ur_ l/2 (t) - (v+ 1/2 (t) - V+ 1/2 (t))) = 0. 

Now, coming back to the conservative variables, we obtain for every i £ [[1, N x ^ and 
every t > 0, 

'w^ x ^ = (vi/2<*) " V1/2W + ^K" + i/ 2 W - ^ r -i/ 2 (*))) 

+ 2A¥ (<-i/ 2 W " <i/ 2 W + ¥>«i/ 2 (*) - <i/ 2 (*))) . 
< u(0, Xi ) = u {xi), ( 5 - 17 ) 
u(t,a) = u a (t), 
u{t,b) = u b {t). 

We recall that by formally setting s = in the scheme (15.511 . we have i> = — ^d x w 
and if = Therefore we can compute 

V tl/2 = -\( d * W t+l/2 atld W tl/2 = 0(«£l/2)' 

where, ^J.i/ 2 > ^ performed using again an ENO reconstruction, see formulae (17.4I) - 
(I7.5I) in Annex ITT1 while the derivatives of w^_ 1 , 2 are approximated using a recon- 
struction polynomial with a centered stencil, see formula (l7.10l) - d7~T2l) in Annex [7721 

We wish to emphasize that the scheme of system ( 15.91 ) reduces to the time advance- 
ment of the single variable u solution of i5. It . 
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5.3. Time discretization 

In order to have a fully discrete scheme, we still need to specify the time discretization. 
According to IfTTIl . we use a discretization based on an explicit Runge-Kutta scheme, 
see 136*1 . for instance. 

We start discretizing the system (15.171) using, for simplicity, a uniform time step At. 
For every i £ [[1, A^]], we denote by ti™ the numerical approximation of u(t m ,Xi) 
with t m = mAt, m = 0, . . . , N t , where N t is the number of time steps. 

The z/-stage explicit Runge-Kutta scheme with v > 1, associated to (15.171 1 can be 
written for every i £ [[1, N x \ as follows, 



u m+i^ u m__J2~ bk F^, (5.18) 
k=i 



where, A = ^ and the stage values are computed at each time step t m and for every 
k epji/J as 



F i = V i+l/2 ~ V i-l/2 + ^K+l/2 " ^-l/ 2 ) " U i-l/2 + V i+l/2 ~ ^1+1/2 ~ W l-l/2 
— U i " ' ' i+l/2 _ 2^ a ^ W /j+1/2' — 



Here (a^, 6^) is a pair of Butcher's tableaux J2TJ, of diagonally explicit Runge-Kutta 
schemes. This finally completes the description of the deterministic numerical method. 

Remark 5.3. In the case when /3 is Lipschitz but possibly degenerate, the authors 
of ifTTl , showed the L 1 -convergence of a semi-discrete in time relaxed scheme, see 
Theorem 1, Section 3 in IfTTl . In fact, they extended the proof of (8) to the case of a v- 
stages Runge-Kutta scheme. Moreover, IfTTIl provided the following stability condition 
of parabolic type, 

At < CAx 2 , (5.20) 

where, C is a constant depending on j3. At the best of our knowledge, no such results 
are available in the case where ft is not Lipschitz. 



(5.19) 



6. Numerical experiments 

We use a Matlab implementation to simulate both the deterministic and probabilistic 
solutions. Concerning the plug-in bandwidth selection procedure described in Section 
|4]and based on J42j, we have improved the code produced by J. S. Matron and available 
on |http : / /www . stat . unc . edu/ f aculty/marron/marron_sof tware . html[ 
by speeding up the root-finding algorithm used to solve (14.71 ) - Furthermore, the deter- 
ministic numerical solutions are performed using the ENO spatial reconstruction of 
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order 3 and a third order explicit Runge-Kutta scheme for time stepping. We point out 
that the deterministic time step, denoted from now on by At^et, is chosen with respect 
to the stability condition ( 15.201 1 . 

6.1. The Classical porous media equation 

We recall that when fl{u) = u.\u\ m ~ l , m > 1, the PDE in (11-11 ) is nothing else but 
the classical porous media equation (PME). The first numerical experiments discussed 
here, will be for the mentioned f3. Indeed, in the case when the initial condition uq is 
a delta Dirac function at zero, we have an exact solution provided in J4), known as the 
density of Barenblatt-Pattle and given by the following explicit formula, 

U(t,x) = t~ fi (C- Kx 2 t~ 2fl )f^, xeR,t>0, (6- 1 ) 

where 



P = — — T' K= ^l — ZT^' C= \ , 7m= / cos X — i 

m+1 2{m+\)m \lm J J-s. 



We would now compare the exact solution (16.11 1 to an approximated probabilistic solu- 
tion. However, up to now, we are not able to perform an efficient bandwidth selection 
procedure in the case when the initial condition is the law of a deterministic random 
variable. Since we are nevertheless interested in exploiting d6.lt . we considered a time 
translation of the exact solution U defined as follows 

v(t,x) = U(t + l,x) Vi6l, Vt>0. (6.2) 



Note that one can immediately deduce from (16.21) . that v still solves the PME but now 
with a smooth initial condition given by 

v (x) = U(l,x) VieR. (6.3) 

In fact, in the case when the exponent m is equal to 3, the exact solution v of the PME 
with initial condition vq(x) = U(l,x) is given by the following explicit formula, 



v(t, x) 



otherwise. 



Simulation experiments: we first compute both the deterministic and probabilistic 
numerical solutions over the time-space grid [0.1.5] x [—2.5,2.5], with space step 
Ax = 0.02. We set At det = 4 x 10" 6 , while, we use n = 50000 particles and a time 
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step At = 2 x 10~ 4 , for the probabilistic simulation. Figures [T](a)-(b)-(c)-(d), display 
the exact and the numerical (deterministic and probabilistic) solutions at times t = 0, 
t = 0.5, t = 1 and t = T = 1.5 respectively. The exact solution of the PME, defined 
in ( 16.41) . is depicted by solid lines. 

Besides, Figure[T](e) describes the time evolution of both the discrete L 2 determin- 
istic and probabilistic errors on the time interval [0, 1.5]. 




0.015 



0.01 



0.005 - 



0.5 



1.5 



Figure 1. - Deterministic (doted line), probabilistic (dashed line) and exact solutions 
(solid line) values at t=0 (a), t=0.5 (b), t=l (c) and t=1.5 (d). The evolution of the L 2 de- 
terministic (dote line) and probabilistic (dashed line) errors over the time interval [0, 1 .5] 
(e). 



The L l errors behave very similarly as well in the present case as in the Heaviside 
case, treated in subsection [O] 
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6.2. The Heaviside case 

The second family of numerical experiments discussed here, concerns /? defined by 
(11.31) . Since we do not have an exact solution of the diffusion problem dl.lt . for the 
mentioned (3, we decided to compare the probabilistic solution to the approximation 
obtained via the deterministic algorithm described in Section|5j Indeed, we shall simu- 
late both solutions according to several types of initial data uo and with different values 
of the critical threshold u c . 

Empirically, after various experiments, it appears that for a fixed threshold u c , the 
numerical solution approaches some limit function which seems to belong to the "at- 
tracting" set 

J = {/ G L l (R)\ J f(x)dx = 1, |/| < u c }; (6.5) 

in fact J is the closure in L 1 of Jo = {/ : R -> K+| (3(f) = 0}. At this point, the 
following theoretical questions arise. 

(1) Does indeed u(t, ■) have a limit Uoo when t — ¥ oo? 

(2) If yes does belong to Jl 

(3) If (2) holds, do we have u(t, ■) = Woo for t larger than a finite time r? 

A similar behavior was observed for different (3 which are strictly increasing after 
some zero. 



6.2.1 Trimodal initial condition 

For the (3 given by (11.31) . we consider an initial condition being a mixture of three 
Gaussian densities with three modes at some distance from each other, i.e. 

u o(x) = 3 {p{x,m,ai) +p(x,fM 2 , 02) +p(x,n 3 ,a 3 )) , (6.6) 

where, 

1 ( x _ n) 2 

p(x, n, a) = — exp( — ). (6.7) 

v27rcr 2<T 

Simulation experiments: for this specific type of initial condition uo, we consider 
two test cases depending on the value taken by the critical threshold u c . We set, for 
instance, fi\ = —fi^ = —4, fi 2 = and o\ = 0.1, 02 = 0.2, 173 = 0.3. 

Test case 1 : we start with u c = 0.15, and a time-space grid [0, 0.6] x [—7, 7], with a 
space step Ax = 0.02. For the deterministic approximation, we set Atjet = 4 x 10 -6 . 
The probabilistic simulation uses n = 50000 particles and a time step At = 2 x 10 . 
Figures|2](a)-(b)-(c), displays both the deterministic and probabilistic numerical solu- 
tions at times t = 0, t = 0.3 and t = T = 0.6, respectively. On the other hand, the 
time evolution of the Z 2 -norm of the difference between the two numerical solutions 
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is depicted in Figure[2](d). 

Test case 2 : we choose now as critical value u c = 0.08 and a time-space grid 
[0,4] x [-8.5,8.5], with a space step Ax = 0.02. We set At det = 4 x 10" 6 and 
the probabilistic approximation is performed using n = 50000 particles and a time 
step At = 2 x 10 -4 . Figures [3](a)-(b)-(c) and[3j(d), show respectively the numerical 
(probabilistic and deterministic) solutions and the L 2 -norm of the difference between 
the two. 

6.2.2 Uniform and Normal densities mixture initial condition 

We proceed with j3 given by dl.31) . We are now interested in an initial condition uo, 
being a mixture of a Normal and an Uniform density, i.e., 

= \ [p{x, -1,0.2) + l [0 ,i]0)) , (6.8) 

where, p is defined in ( 16.71 . 
Simulation experiments: 

Test case 3 : we perform both the approximated deterministic and probabilistic 
solutions in the case where u c = 0.3, on the time-space grid [0, 0.5] x [—2.5, 2], with 
a space step Ax = 0.02. We use n = 50000 particles and a time step At = 2 x 10~ 4 , 
for the probabilistic simulation. Moreover, we set At^et = 4 x 10~ 6 . Figures [4] (a) - 
(b)-(c) illustrate those approximated solutions at times t = 0, t = 0.1 and t = T = 
0.5. Furthermore, we compute the L 2 -norm of the difference between the numerical 
deterministic solution and the probabilistic one. Values of this error, are displayed in 
Figure |4](d), at each probabilistic time step. 

6.2.3 Uniform densities mixture initial condition 

Now, with f3 given by (11.31 1. we consider an initial condition uq being a mixture of 
Uniform densities, i.e., 

1 3 5 

u (x) = -l[ 0)1 ]O) + -l [ _i ) i ] (ar) + -t { 6 i2] (x), (6.9) 

Simulation experiments: 

Test case 4 : we approximate the deterministic and probabilistic solutions in the 
case where u c = 0.3, on the time-space grid [0,0.6] x [—1.5,3.5], with a space step 
Ax = 0.02. The deterministic time step Atdet = 4 x 10~ 6 , while the probabilistic 
solution is computed using n = 50000 particles and a time step At = 2 x 10 -4 . 
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We illustrate in Figure[5](a)-(b)-(c), both the deterministic and probabilistic numerical 
solutions at times t = 0, t — 0.1 and t = T = 0.6 respectively; while the time 
evolution of the L 2 -norm of the difference between them, is shown in Figure[5](d). 

6.2.4 Square root initial condition 

Finally, the last test case concerns an initial condition uq defined as follows : 



Simulation experiments: 

Test case 5 : we simulate the probabilistic and deterministic solutions over the time- 
space grid [0,0.45] x [—2,2], using a space step Ax = 0.02 and setting the critical 
threshold u c = 0.35. Moreover, the deterministic time step A^ ei = 4 x 10~ 6 . On the 
other hand, we use n = 50000 particles and a time step At = 2 x 10~ 4 to compute the 
probabilistic approximation. 

Figures |6](a)-(b)-(c), show both the deterministic and probabilistic numerical solu- 
tions at times t = 0, t = 0.04 and t — T — 0.45, respectively. 

The evolution of the L 2 -norm of the difference between these two solutions, over 
the time interval [0,0.45], is depicted in Figure[6](d). 

6.3. Concluding remarks 

(1) Figure|7l displays a single trajectory for each one of the test cases described above. 

In fact, we observe that, in all cases, the process trajectory stops not later than 
the instant of stabilization of the macroscopic distribution. 

(2) We have performed deterministic and probabilistic numerical solutions for dl.lt . 

with a coefficient fj defined by ( 11.3b . Even though the procedures being used 
were different, the simulation experiments clearly show that both methods pro- 
duce very close approximated solutions, all over the considered time interval. 

(3) We point out that, the error committed by the Monte Carlo simulations largely 

dominates the one related to the Euler scheme. Consequently, the choice of the 
probabilistic time step is not so important. 

(4) The probabilistic algorithm can be parallelized on a Graphical Processor Unit 

(GPU), such that we can speed-up its time machine execution; on the other 
hand, for the deterministic algorithm, this transformation is far from being ob- 
vious, see |[T8l . 




(6.10) 
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Figure 2. - Test case 1: Deterministic (solid line) and probabilistic (doted line) solution values at t=0 
(a), t=0.3 (b), t=0.6 (c). The evolution of the L"-norm of the difference over the time interval [0, 0.6] (d). 
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Figure 3. - Test case 2: Deterministic (solid line) and probabilistic (doted line) solution values at t=0 
(a), t=2 (b), t=4 (c). The evolution of the L 2 -norm of the difference over the time interval [0, 4] (d). 
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Figure 5. - Test case 4: Deterministic (solid line) and probabilistic solution (doted line) values at t=0 
(a), t=0.1 (b), t=0.6 (c). The evolution of the L 2 -norm of the difference over the time interval [0, 0.6] (d). 
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Figure 7. - Representation of a process trajectory for the Test case 1 (a), Test case 2 (b), Test case 3 (c), 
Test case 4 (d) and Test case 5 (e), respectively. 
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7. Annexes 

Let V G R Nx such that Vi = v(xi), Vi G [[1,N X }\ where v is a function defined on 
[a, b}. Note that the points (xi) are still defined as in Section[5] Moreover, A4 m ,n(J&) 
denotes the linear space of real matrices with m rows and n columns. 

7.1. Interpolating polynomial of a function 

We aim to approximate v(x i+1 / 2 ) an d K^i-l 12) f° r every i G [1, N x \ ■ In order to do 
this, we use properly chosen Lagrange interpolation polynomials of degree k — 1 . 

On every interval (or cell) Ii = [2^-1/2, x i+ i/ 2 ], with i G p, A^J , we construct 
an interpolation polynomial , by selecting fe consecutive points containing X{ : the 
so-called stencil denoted by 

^ W — {-^i— rj • • • ) -^i+s} 

and defined by {xi_ r , Xi_ r+ i , . . . , Xj+ s -i , Xi +S }, where r, s are positive integers and 
r + s + 1 = k. We denote by the value taken by r for the interval Ii with an 
ENO stencil, see (44). 

The Lagrange interpolation polynomial of degree k— 1, on the interval associated 
to the stencil S(i) is then given by : 

fe-i 

p S( a: )=Z) yi -^ L i 1 ( !r )' Va;e/i ' (7 - 1} 

j=o 

where, 

4 ] (ar) = J] X ~ Xl J r+l . (7.2) 
¥i 

Now, we need to compute the polynomial defined in d7.2l l at the points x^i/2 and 
Xi+i/2- In fact, since the points are equidistant, we have for every (r,j) G [[0, k — l]f , 

^■ (^-1/2) = n — a _ r ' 

¥i 

< 

^■(^+1/2) = n — a _ / • 

^ ¥i 
Then, we define C G fc(R), as follows 

fc— l 

C r+ i, i+1 = II r ~^r/ /2 ' V(r,j)G[[0,A:]]x[[0,fc-l]]. (7.3) 
¥i 
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Substituting d73} in (TTTT) and using the ENO stencil, we get Mi G [1, N x \ 

fc-l 

v( Xi _ l/2 ) » w + _ 1/2 = Pg (x^t/a) = Fi-fl^+jC^+jj+i, (7.4) 

i=o 
fe-i 

v{x i+l/2 ) » «^ 1/2 = p£ ( (z i+ i/ 2 ) = X] ^-RW+j C flW+2J+i- ( 7 - 5 ) 

7.2. Interpolation polynomial for the derivative of a function 

Now, we would like to approximate ^(x j), ^7(2^-1/2) and ^(2^+1/2), f° r every 
i G [[1, iVa,]] . In fact, deriving equation (17.11) . implies 



(x) = 5^y i _ r+j — i-(x), Vxg/,. (7.6) 



dx ^— ' 

On the other hand, for every j G [[0, / — 1]] , we have 

fc-i fc-i 

XI ] [ ( X ~ X i~r+l) 



jt-M m=0 1=0 

"^-W = fcZl . G (7.7) 

^J(^»— T+j %i— r+l) 
1=0 

Hi 



Since the points are equidistant, we get 

fc-l fc-l 

e n ( r -o 



dL 



r] m =0 1=0 

m^j Hj,m 



H*0 



and 



dx v *' fc-l 

AxH(j-l) 

1=0 

¥i 



fc-i fc-i fc-i fc-i 

£ (r- 1-1/2) £ (r- 1 + 1/2) 

m=0 Z=0 jr[r] m=0 1=0 



jrl »J m=u t=u ///"-- 



■(Si-1/2) = TZ1 . -7-(^+l/2) 



fc-l ' dx v-n-V^ fc _j 



z=o z=o 
¥j ¥i 
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Then, we define B e 



fc-i fc-i 

e nc r -o 

m=0 1=0 

Wi = !!!f n2f — > v ^') e I - fc - l ti > (7 - 8) 



/-o 



and D £ .Mfc+i,fc(R) as follows 
fc-i fe-i 

E U<r-i-V2) 

m=0 1=0 

B r+lj+1 = m ^ ¥j '7_ 1 , V(r,j)e[[0,fc]]x[[0,fc-1]]. (7.9) 

Z=0 

¥i 

Therefore, replacing (17.81) and (17.91 1 in (17.61) . for every i 6 [[1 , iV x ]] , we obtain : 

— (Xi) W dUj = — ~(Xi) = Vi-r+jDr+ij + i, (7.10) 

x x i=o 

^(^-1/2) ~ dU z-l/2 = -^(^-1/2) = E^- r +J" ID)r + 1 'J +1 ' (? - U) 
X X 3=0 

dv cflP^- k ~ l 

^(^+1/2) ~ d«- 1/2 = -^(i i+1/2 ) = 5^Vi_ r+i D r+ . 2)J - + i. (7.12) 



A probabilistic algorithm approximating a singular PDE 



41 



Acknowledgments. Part of the work was done during the stay of the first and third 
named authors at the Bielefeld University (SFB 701 and BiBoS). They are grateful for 
the invitation. 



References 

1. D. Aregba-Driollet, R. Natalini and S. Tang, Explicit diffusive kinetic schemes for nonlin- 
ear degenerate parabolic systems, Math. Comp. 73 (2004), pp. 63-94 (electronic). 

2. P. Bak, How Nature Works: The science of Self-Organized Criticality. Springer- Verlag 
New York, Inc, 1986. 

3. V. Barbu, M. Rockner and F. Russo, Probabilistic representation for solutions of an irreg- 
ular porous media type equation: the irregular degenerate case. To appear: Prob. Th. Rel. 
Fields. Available at |http : //hal . inria . f r/inria-00410248/f r71 

4. G. I. Barenblatt, On some unsteady motions of a liquid and gas in a porous medium, Akad. 
Nauk SSSR. Prikl. Mat. Meh. 16 (1952), pp. 67-78. 

5. S. Benachour, P. Chassaing, B. Roynette and P. Vallois, Processus associes a Vequation 
des milieux poreux, Ann. Scuola Norm. Sup. Pisa CI. Sci. (4) 23 (1996), pp. 793-832 
(1997). 

6. P. Benilan, H. Brezis and M. G. Crandall, A semilinear equation in L l (M. N ), Ann. Scuola 
Norm. Sup. Pisa CI. Sci. (4) 2 (1975), pp. 523-555. 

7. P. Benilan and M. G. Crandall, The continuous dependence on <p of solutions of u t — 
Aip(u) = 0, Indiana Univ. Math. J. 30 (1981), pp. 161-177. 

8. A. E. Berger, H. Brezis and J. C. W. Rogers, A numerical method for solving the problem 
u t - A/0) = 0, RAIRO Anal. Numer. 13 (1979), pp. 297-312. 

9. P. Blanchard, M. Rockner and F. Russo, Probabilistic representation for solutions of an 
irregular porous media type equation, Ann. Probab. 38 (2010), pp. 1870-1900. 

10. M. Bossy and D. Talay, A stochastic particle method for some one-dimensional nonlinear 
p.d.e, Math. Comput. Simulation 38 (1995), pp. 43-50. Probabilites numeriques (Paris, 
1992). 

11. , A stochastic particle method for the McKean-Vlasov and the Burgers equation, 

Math. Comp. 66 (1997), pp. 157-192. 

12. F. Bouchut, F. R. Guarguaglini and R. Natalini, Diffusive BGK approximations for nonlin- 
ear multidimensional parabolic equations, Indiana Univ. Math. J. 49 (2000), pp. 723-749. 

13. A. W. Bowman, An alternative method of cross-validation for the smoothing of density 
estimates, Biometrika 71 (1984), pp. 353-360. 

14. H. Brezis and M. G. Crandall, Uniqueness of solutions of the initial-value problem for 
u t - A(p(u) = 0, J. Math. Pures Appl. (9) 58 (1979), pp. 153-163. 

15. R. Cafiero, V. Loreto, L. Pietronero, A. Vespignani and S. Zapperi, Local regidity and 
self-organized criticality for avalanches, Europhysics Letters 29 (1995), pp. 1 1 1-1 16. 

16. P. Calderoni and M. Pulvirenti, Propagation of chaos for Burgers' equation, Ann. Inst. H. 
Poincare Sect. A (N.S.) 39 (1983), pp. 85-97. 



42 



Nadia Belaribi, Francois Cuvelier and Francesco Russo 



17. F. Cavalli, G. Naldi, G. Puppo and M. Semplice, High-order relaxation schemes for non- 
linear degenerate diffusion problems, SIAM J. Numer. Anal. 45 (2007), pp. 2098-2119 
(electronic). 

18. F. Cuvelier, Implementing Kernel Density Estimation on GPU: application to a proba- 
bilistic algorithm for PDEs of porous media type, Technical report. In preparation. 

19. J. M. Dawson, Particle simulation of plasmas, Rev. Modern Phys. 55 (1983), pp. 403-447. 

20. A. Figalli and R. Philipowski, Convergence to the viscous porous medium equation and 
propagation of chaos, ALEA Lat. Am. J. Probab. Math. Stat. 4 (2008), pp. 185-203. 

21. E. Hairer, S. P. N0rsett and G. Wanner, Solving ordinary differential equations. I, sec- 
ond ed, Springer Series in Computational Mathematics 8. Springer- Verlag, Berlin, 1993, 
Nonstiff problems. 

22. A. Harten and S. Osher, Uniformly high-order accurate nonoscillatory schemes. I, SIAM 
J. Numer. Anal. 24 (1987), pp. 279-309. 

23. R. W. Hockney and J. W. Eastwood, Computer simulation using particles. McGraw-Hill, 
New York, 1981. 

24. S. Jin and C. D. Levermore, Numerical schemes for hyperbolic conservation laws with 
stiff relaxation terms, J. Comput. Phys. 126 (1996), pp. 449-467. 

25. S. Jin and Z. P. Xin, The relaxation schemes for systems of conservation laws in arbitrary 
space dimensions, Comm. Pure Appl. Math. 48 (1995), pp. 235-276. 

26. M. C. Jones, J. S. Marron and S. J. Sheather, A brief survey of bandwidth selection for 
density estimation, J. Amer. Statist. Assoc. 91 (1996), pp. 401-407. 

27. B. Jourdain, Probabilistic approximation for a porous medium equation, Stochastic Pro- 
cess. Appl. 89 (2000), pp. 81-99. 

28. B. Jourdain and S. Meleard, Propagation of chaos and fluctuations for a moderate model 
with smooth initial data, Ann. Inst. H. Poincare Probab. Statist. 34 (1998), pp. 727-766. 

29. J. Kacur, A. Handlovicova and M. Kacurova, Solution of nonlinear diffusion problems by 
linear approximation schemes, SIAM J. Numer. Anal. 30 (1993), pp. 1703-1722. 

30. I. Karatzas and S. E. Shreve, Brownian motion and stochastic calculus, second ed, Grad- 
uate Texts in Mathematics 113. Springer- Verlag, New York, 1991. 

31. H. P. Jr. McKean, Propagation of chaos for a class of non-linear parabolic equations., 
Stochastic Differential Equations (Lecture Series in Differential Equations, Session 7, 
Catholic Univ., 1967), Air Force Office Sci. Res., Arlington, Va., 1967, pp. 41-57. 

32. S. Meleard and S. Roelly-Coppoletta, A propagation of chaos result for a system of parti- 
cles with moderate interaction, Stochastic Process. Appl. 26 (1987), pp. 317-332. 

33. K. Oelschlager, A law of large numbers for moderately interacting diffusion processes, Z. 
Wahrsch. Verw. Gebiete 69 (1985), pp. 279-322. 

34. , A fluctuation theorem for moderately interacting diffusion processes, Probab. 

Theory Related Fields 74 (1987), pp. 591-616. 

35. , Simulation of the solution of a viscous porous medium equation by a particle 

method, SIAM J. Numer. Anal. 40 (2002), pp. 1716-1762 (electronic). 

36. L. Pareschi and G. Russo, Implicit-Explicit Runge-Kutta schemes and applications to hy- 
perbolic systems with relaxation, J. Sci. Comput. 25 (2005), pp. 129-155. 

37. E. Parzen, On estimation of a probability density function and mode, Ann. Math. Statist. 
33 (1962), pp. 1065-1076. 



A probabilistic algorithm approximating a singular PDE 



43 



38. R. Philipowski, Interacting diffusions approximating the porous medium equation and 
propagation of chaos, Stochastic Process. Appl. 1 17 (2007), pp. 526-538. 

39. I. S. Pop and W. Yong, A numerical approach to degenerate parabolic equations, Numer. 
Math. 92 (2002), pp. 357-381. 

40. M. Rudemo, Empirical choice of histograms and kernel density estimators, Scand. J. 
Statist. 9 (1982), pp. 65-78. 

41. D. W. Scott and G. R. Terrell, Biased and unbiased cross-validation in density estimation, 
J. Amer. Statist. Assoc. 82 (1987), pp. 1131-1146. 

42. S. J. Sheather and M. C. Jones, A reliable data-based bandwidth selection method for 
kernel density estimation, J. Roy. Statist. Soc. Ser. B 53 (1991), pp. 683-690. 

43. R. E. Showalter, Monotone operators in Banach space and nonlinear partial differential 
equations, Mathematical Surveys and Monographs 49. American Mathematical Society, 
Providence, RI, 1997. 

44. C. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for 
hyperbolic conservation laws, Advanced numerical approximation of nonlinear hyper- 
bolic equations (Cetraro, 1997), Lecture Notes in Math. 1697, Springer, Berlin, 1998, 
pp. 325^132. 

45. B. W. Silverman, Density estimation for statistics and data analysis, Monographs on 
Statistics and Applied Probability. Chapman & Hall, London, 1986. 

46. D. W. Stroock and S. R. S. Varadhan, Multidimensional diffusion processes, Classics in 
Mathematics. Springer-Verlag, Berlin, 2006, Reprint of the 1997 edition. 

47. A. S. Sznitman, Topics in propagation of chaos, Ecole d'Ete de Probabilites de Saint-Flour 
XIX— 1989, Lecture Notes in Math. 1464, Springer, Berlin, 1991, pp. 165-251. 

48. G. R. Terrell, The maximal smoothing principle in density estimation, J. Amer. Statist. 
Assoc. 85 (1990), pp. 470-477. 

49. G. R. Terrell and D. W. Scott, Oversmoothed nonparametric density estimates, J. Amer. 
Statist. Assoc. 80 (1985), pp. 209-214. 

50. M. P. Wand and M. C. Jones, Kernel smoothing, Monographs on Statistics and Applied 
Probability 60. Chapman and Hall Ltd., London, 1995. 

51. M. Woodroofe, On choosing a delta-sequence, Ann. Math. Statist. 41 (1970), pp. 1665- 
1671. 

52. J. F. Grotowski Z. I. Botev and D. P. Kroese, Kernel density estimation via diffusion, 
Submitted to the Annals of statistics. (2007). 



Author information 

Nadia Belaribi, Laboratoire d'Analyse, Geometrie et Applications (LAGA), Universite Paris 
13, 99, avenue Jean-Baptiste Clement, F-93430 Villetaneuse and ENSTA ParisTech, Unite de 
Mathematiques appliquees, 32, Boulevard Victor, F-75739 Paris Cedex 15, France. 
Email: belaribi@math . univ-parisl3 . f r 

Francois Cuvelier, Laboratoire d'Analyse, Geometrie et Applications (LAGA), Universite 
Paris 13, 99, avenue Jean-Baptiste Clement F-93430 Villetaneuse, France. 
Email: cuvelier@math . univ-parisl3 . f r 



44 



Nadia Belaribi, Francois Cuvelier and Francesco Russo 



Francesco Russo, ENSTA ParisTech, Unite de Mathematiques appliquees, 32, Boulevai'd Vic- 
tor, F-75739 Paris Cedex 15, INRIA Rocquencourt and Cermics Ecole des Ponts et Chaussees, 
Projet MATHFI Domaine de Voluceau, BP 105 F-78153 Le Chesnay Cedex, France. 
Email: francesco.russo@ensta-paristech.fr 



